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Abstract 



o 

C^ ' Channelrhodopsins-2 (ChR2) are a class of light sensitive proteins that offer the ability to use light stimulation 

O |. to regulate neural activity with millisecond precision. In order to address the limitations in the efficacy of the 

wild-type ChR2 (ChRwt) to achieve this objective, new variants of ChR2 that exhibit fast mono-exponential 

o: 

Cn| . photocurrent decay characteristics have been recently developed and validated. In this paper, we investigate 

whether the framework of transition rate model with 4 states, primarily developed to mimic the bi-exponential 

^T' . photocurrent decay kinetics of ChRwt, as opposed to the low complexity 3 state model, is warranted to 

mimic the mono-exponential photocurrent decay kinetics of the newly developed fast ChR2 variants: ChETA 



(Gunaydin et al., Nature Neurosci, 13:387-392, 2010) and ChRET/TC (Bcrndt et al., PNAS, 108:7595-7600, 
2011). We begin by estimating the parameters for the 3-state and 4-state models from experimental data 
on the photocurrent kinetics of ChRwt, ChETA and ChRET/TC. We then incorporate these models into 
*/~) ' a fast-spiking intcrncuron model (Wang and Buzsaki., J Neurosci, 16:6402-6413.1996) and a hippocampal 

m ; 

^sO , pyramidal cell model (Golomb et al., J Neurophysiol, 96:1912-1926, 2006) and investigate the extent to 

in 



which the experimentally observed neural response to various optostimulation protocols can be captured by 
these models. We demonstrate that for all ChR2 variants investigated, the 4 state model implementation 
is better able to capture neural response consistent with experiments across wide range of optostimulation 
protocol. We conclude by analytically investigating the conditions under which the characteristic specific 
to the 3-state model, namely the mono-exponential photocurrent decay of the newly developed variants of 
ChR2, can occurs in the framework of the 4-state model. 
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1. Introduction 



Optogenetics, an emerging optical neurostimulation technique, employs light activate d ion channels to 



excit e or suppress impulse activity in neurons with high temporal and spatial resolution (jDcisseroth et al 



20061 ) . The use of light act ivation for modulating properties of neurons, an idea that goes back to Nobelist 



Francis Crick (Crick 



1979 ). has generated great excitement in the biological science community ([Kuehn 



2010J). Indeed, in recent years there has been a great s urge in research efforts using this technology to 



addre ss fundamental questions in the field of neuroscience ( Huber et al 



20101 ) and neurological disease (jGradinaru et al 



2009; 



T0nnesen et al. 



2008 



2009; 



Cardin et al 



Kravitz et al 



2009: 



2010) 



Bass et al 



The basic idea behind optogcnetic technology is to enable the induction of light activated ion channel 
proteins (opsins) into neuronal membranes of living animals using techniques from molecular and genetic 
engineering. The activity of the opsin containing neurons can then be controlled directly via light stimu- 
lation. Specifically, experimental strategies have been developed for delivering channclrhodopsin-2 (ChR2), 
a subfamily of blue light gated-cation opsins, and halorhodopsin (Halo/NpHR), a subfamily of yellow light 
gatcd-anion opsins, such th at optostimulation can either activate (generate action potent ial) or deactivate 



(hyperpolarize) the neuron (jBovden et al 



2005 



Deisseroth et al 



2006 



Zhang et al 



20071) . A key merit of 



this technique is that it offers the ability to regulate neuronal activity with millisecond precision, which in 
turn allows for fine control of neuronal activity patterns in the brain region of interest. In additio n, these 



20091) . Thus, 



proteins can be engineered to be expressed only in certain types of neurons (| Cardin et al.l 
this technique offers capability to control neuronal activity with high degree of temporal accuracy in a 
cell specific manner, a significant advantage over traditional techniques such as electrical stimulation and 
pharmacological approaches. 

Recent in vitro experiments to characterize the response of neurons expressing the wild-type ChR2 
(ChRwt) to various optostimulation protocols have highlig hted several limitation s in the precision of op- 



20101) . Specifically, short light 



togenetic control that can be achieved via optostimulation (IGunavdin et al. 
pulse (2 ms) stimulation of neurons expressing ChRwt can result in the generation of extra spikes; fast- 
spiking intcrncurons driven with periodic short light pulse stimulation in the gamma range (30-80 Hz) fail to 
respond with spiking in the gamma range while higher frequency optostimulation (>80 Hz) triggers plateau 
potential neural response. In order to overcome these limitations, ChR2 mutations with fast photocurrcnt 
decay ch arcteristics have been engineered and vali dated, which offer enhanced and precise control of neural 



activity (IGunavdin et al 



2010: 



Berndt et al 



2011) 



From a modeling perspective, two transition rate models, a 3-state (JNikolic et al 



20061 ) and a 4-state 



model (jNikolic et al.ll2009J ). are currently available to mimic the ChRwt photocurrent kinetics. The 4 state 
model was proposed to specifically capture the bi-expon ential decay profile of the ChRwt photocurrent 



20091 ). However, to our knowledge, 



following the termination of the light stimulation protocol (jNikolic et al. 
the bi-exponential photocurrent decay characteristic has not been reported in other studies using ChR2wt or 
any other ChR2 variant for optogenetic manipulation of neural activity. Therefore, the benefit of enhanced 
repertoire of dynamical behavior offered by an additional state variable in the transition rate model relative 
to the increase in the complexity of the model is called into question. Moreover, it is unclear whether 
either of these transition rate models can capture th e photocurrent kinetics of the recently engineered ChR2 



variants (jGunavdin et al. 



2010: 



Berndtetal 



20111) . Furthermore, no study has yet investigated whether 
neural activity elicited in model neurons using photocurrents generated by these transition rate models is in 
agreement with experimental results obtained when various optostimulation protocols are delivered. Given 
that modern control engineering analysis and design tools depend crucially on mathematical models of the 
actuation system (ChR2 photocurrents) and the phenomena to be controlled (neural network activity), it is 
of vital importance to characterize the efficacy of mathematical models of ChR2 photocurrents to reproduce 
the main features of neural response to various optostimulation protocols. 

This study is focused on a systematic investigation of the ability for a 3-state and a 4-state transition 
rate model respectively, to mimic photocurrent characteristics of ChRwt that exhi bits mono-exponentia l 



photocurrent dec ay and the two new ly engineered ChR2 variants, namely, ChETA (IGunavdin et al 



and ChRET/TC (JBerndt et al 



2010) 



20111 ) in order to determine which of the models can effectively capture the 



the dynamical response of the neurons expressing these variants to various optostimulation protocols. We 
begin by first identifying the parameters for the 3-state and the 4-state model of each of the three ChR2 
variants. We then incorporate these mo dels into two well e stablished neuron models, namely the Golomb 



model for hippocampal py ramidal cells (Golomb et al 



hippocampal interneurons (jWang and Buzski 



20061 ) and the Wang-Buszaki model for fast spiking 



19961) and study the dynamical response of these neuron models 
to different optostimulation protocols. Finally, we present analytical results on the conditions under which 
a characteristic specific to the 3-state model, namely the mono-exponential decay of the ChR2 photocurrent 
kinetics, can occur in the 4-state model. 

The paper is organized as follows: We first summarize the key experimental results incorporated in 
our development of the 3-state and 4-state transition rate models of ChR2 photocurrent kinetics. We then 
present the mathematical framework of the transition rate models, which is then followed by the presentation 
of the neuron models studied in this manuscript. Results from our analysis of the two transition rate models 



Experimental Data 


ChR2 Variant 


(ms) 


T rise (**) 

(ms) 


Tin 

(ms) 


Toff 

(ms) 


(ms) 


R 


-'peak 

(nA) 


Vhold 

(mV) 


flight 
/ mW \ 


Gunaydin et al. 
(Gunavdin et al. 2010) 


ChRwt 
ChETA 


2.4 
0.9 


0.2 
0.08 


55.5 
15 


9.8 
5.2 


10700 
1000 


0.4 
0.6 


0.848 
0.645 


-100 


50 


Berndt et al. 
(Berndt et al. 2011) 


ChRwt 
ChRET/TC 


2.65 
2.17 


0.23 
0.18 


9.6(*) 
11(*) 


11.1 
8.1 


10700 
2600 


0.27(*) 
0.31(*) 


0.967 
1.420 


-75 


42 



Table 1: Experimental Data. The values of various time constants and photocurrent quantifiers determined experimentally 
for ChRwt and the two fast variants ChETA and ChRET/TC are presented. These data have been employed in the evaluation 
process of the parameters for the 3 and 4 state ChR2 model for each variant. (*) These data have been obtained by direct 
correspondence with the authors of the experimental study. (**) These data have been analytically evaluated (see Appendix, 
Section 6.4.1) 



are then presented beginning with results on our estimates of parameter values for the two transition rate 
models. The discussion section then follows, where we explain the relevance of this study for the development 
of neural control strategies using optostimulation based actuation systems. 

2. Methods 



2.1. Experimental data 

The 3- and 4-state transition rate models for each of the ChR2 variant inve stigated in this 



are designed to m atch the experimental data currently available in the literature IjGunavdin et al 



Berndt et al. 



paper 



2010: 



20111 ). The available data correspond to photocurrent time profiles obtained from whole cell 



patch clamp recordings from neurons expressing ChR2 in response to prolonged (1 s) optostimulation. 

The following empirically estimated parameters derived from the time profile of the recorded ChR2 
photocurrcnts were employed in the development of the two transition rate models: the time t p , of peak 
ChR2 photocurrent (/peak) after the optostimulation pulse was turned on, the time constant T; n , of the 
exponential decay of the ChR2 photocurrent from peak to steady state plateau (ipiat), the time constant 
t s , of the exponential decay of the photocurrent from plateau to zero when light pulse is turned off, the ratio 
R = -/plat //peak and the value of I poa k. The first parameter t p , allows us to estimate the time constant r r i sc , 
of the photocurrent rise (see Appendix, Section 6.4.1 for more details). In addition, the ChR2 photocurrent 
recovery curve (corresponding to the recovery of the peak current following stimulation with a second light 
pulse) provides the recovery time constant T r . The values for these parameters derived from experimental 
data arc summarized in Table [TJ 



2.2. Transition rate models 

ChR2 ph otocurrent k inetics can be described using the general mathematical framework of transition 



rate models (Koch 



19981) . Such models are described in terms of states S\ ■ ■ ■ S n , n > 2, such that the 



transitions between any two states is given as: 



J, T O , 



(l) 



where the rates rij and Tji govern the transition between states Si and Sj . The fraction of channels in the 
state Sj is denoted by Sj and obeys the following first order differential equation 



ab i \r~^ \r~^ 

3=1 J=l 



(2) 



such that X)"_j Si = 1. 

In Figure [TJ we present a schematic diagram of the 3-state and the 4-state transition rate models fo r 



ChR2 photocurrent kinetics based on the framework described above (JNagel et al 



2003: 



Nikolic et al 



200G) 



B 



X 




^J> ^)"^JU 



Closed State 



¥ Open State 



<^ 



• Closed States 
> Open States 



Figure 1: Schematic diagram of the 3 and 4-state models. A. Representation of the 3-state model with an open state 
O, a closed C and a dark state D. B. The 4-state model is represented by transitions among two open states (Oi and O2) and 
two closed states (C\ and C2). 



2.2.1. 3-state transition rate model for ChR2 photocurrent kinetics 



The 3-state model ( Nagel et al. 



2003: 



Nikolic et al 



20061) is schematically described in Figure [T]/L The 



model simulates the photocurrent kinetics of ChR2 through a set of transitions between three distinct states: 
closed C, open O and a desensitized state D. In the absence of optostimulation, ChR2 molecules are assumed 
to be in state C. Upon illumination with light of appropriate wavelength (ss 475 nm), the ChR2 molecules 



undergo conformational changes and transition to state O, which then spontaneously decays into a closed 
but desensitized state D. ChR2 molecules in state D arc not available to photoswitch on optostimulation. 
Finally, following a prolonged recovery time period, which is much slower than the time scales involved in 
the light induced C — > O and the spontaneous O — » D transitions, the protein returns to the conformation 
of the closed state C . 

The light mediated transition process from the closed state C to the open state O is very fast (< 1 
ms). The transition from the open state O t o the desensitized s tate D depends on the pH of the media 



(extracellular solution) and lasts wlO-400 ms (INikolic et al 



20061 ). In our modeling framework we do not 



consider the pH dependence of the transitions; instead, we will rely on the data reported in Table 1, which 
has been collected for a media of pH = 7.5. If c, o and d denote the fraction of ChR2 molecules in each of 
the three states at any given instant in time, then following from equation [21 the transition rate model for 
ChR2 kinetics can be described via the following set of ordinary differential equations: 



6 = P(l - o - d) - G d o 

d = GdO — G r d 



(3) 



where, o + d + c = 1. The parameters: P, the light dependent excitation rate for the C — >• O transition, 
Gd, the rate constant for the O —> D transition and G r , the rate constant for the D — >• C transition, are 
the model parameters to be determined from experimental data summarized in Section 2.1. The ChR2 
photocurrent entering the neuron membrane is given as: 



^ChR2 = giVo 



(4) 



where V is the membrane potential of the neuron expressing ChR2 and gi is the maximal conductance of 
the ChR2 ion channel. 

Closed form analytical expression for o and therefore JchR2 (under voltage clamp conditions, following 
from equation ??) under both c onstant light intensity optostimulation and no light conditions can be obtained 



see (Nikolic et al 



200G, 



20091) for a complete derivation) from equation [3] as: 



°Sau(*) = Cie- Mt + C 2 e-^ + o plat 



(5) 



Experimental Data 


ChR2 Variant 


P (ms- 1 ) 


G d {ms-') 


G r {ms l ) 


gi(*)(IIC) 


gi(*)(SIC) 


Gunavdin (Gunavdin et al. 2010) 


ChRwt 


0.0179 


0.1020 


9.3458e-005 


0.07 


3.687 




ChETA 


0.0651 


0.1923 


le-003 


0.03314 


0.7588 


Berndt (Berndt et al. 2011) 


ChRwt 


0.1048 


0.0901 


9.3458e-005 


0.03256 


3.3728 


ChRET/TC 


0.0895 


0.1235 


3.8462e-004 


0.06097 


1.899 



Table 2: Parameters of the 3-state model. The parameters employed in the 3-state model have been determined for ChRwt 
and the two fast variants ChETA and ChRET/TC using experimental data presented in Table [l] and the formulas ??-?? from 
this section. (*)This parameter is obtained by solving the optimization problem of the photocurrent peak by minimizing the 
error function provided by equation ?? for different initial conditions (IIC = Ideal Initial Conditions; SIC = Special Initial 
Conditions). 



G d t. 



"«(*) = Ce~^ 



(0) 



where the coefficients C\,Ci. the steady state fraction of ChR2 molecules in open state, o v \ a t and the time 
decays rate constants Ai,A2(> Ai) depend on the model parameters {P, Gd,G r } and the initial conditions 
{o U; cio} (see Section 6.1 of the Appendix for the formulas). The transition rate parameters {Gd,G r } and 
the photocurrent decay rate constants {Ai,A2} are directly obtained from the experimental data on ChR2 
photocurrents as: 



Gd — ; G r — — ; Ai — — A2 — 

1~off <r ^"in Trise 



(7) 



In order to estimate P, we use the analytical solution to equation [3] under the light on conditions to obtain 
the closed form equation (see Appendix, Sections 6.1 and 6.2): 



P = Xi 



G r Gd 



Ai — G r — Gd 



(8) 



In Table [21 we summarize the values obtained for the 3-state model parameters for all ChR2 variants from 
equation [7j equation [8] and the corresponding experimental data. 



2.3. J^-state transition rate model for ChR2 photocurrent kinetics 



20091) . This model 



The 4-state transition rate model is schematically described in FigureQjB (JNikolic et al. 
was primarily developed to account for the key experimental findings related to the bi-exponcntial decay 
rates for the wild-type ChR2 ph otocurrent following the termination of p rolonged optostimulation pulse, as 



observed in the experiments of (Ish izuka et al 



2006: 



Nikolic et al 



2009J). The 4-state model simulates the 



photocurrent kinetics for the light activated ChR2 channel through two sets of intratransitional states: Cl-O- 



01, which is more dark adapted and C2 -f-> 02, which is more light adapted. In absence of optostimulation, 
the ChR2 molecules are assumed to be in the closed state CI. In the light adapted state, however, ChR2 
molecules are distributed across all four states, with increasing preference to be in state 02 as the duration of 
optical illumination increases. For a given level of optical excitation, an equilibrium is established between the 
states 01 and CI. As the duration of illumination increases, the transition to state 02 becomes significant. 
The 02 state is also in equilibrium with the state C2, which in turn slowly transition back to the state CI 
following the termination of the optical signal with time scale of the recovery period. 

If the fraction of ChR2 molecules in each of the four states CI, C2, 01 and 02 are given by ci, C2, o\ 
and 02 respectively, then the four state model for the ChR2 photocurrent kinetics can be represented via the 
following set of ordinary differential equations: 



o'i = P\s{\ - c 2 - o\ - o 2 ) - (Gdi + ei 2 )oi + e 2 i.o 2 

oi = P2SC2 + e 12 o 1 - (Gd2 + e2i)o 2 

C 2 = G d2 02 - (P 2 S + G r )c 2 

s = (S (6) - s)/rchR2 



(9) 



such that ci + C2 + 01 + o 2 = 1. The parameters P\ and Pi are the maximum excitation rates of the first 
and second closed states, Gd\ and Gd2 are the rate constants for the OX — > CI and the 02 — > C2 transitions 
respectively, ei 2 and e 2 i are the rate constants for Ol — > 02 and 02 — > Ol transitions respectively and 
G r is the recovery rate of the first closed state after the light pulse is turned off. We note that, while the 
photon ab sorption and t he isomerization of the r e tinal compound in the rhodopsin is almost instantaneous 



(~ 200 fs 



Wang et al 



()l994h : 



Schoenlein et al 



(|l99ll) V the conformational change leading to the open 



state configuration is a slower process. The function s is designed to capture the temp oral kinetics of this 



conformational change in the protein. This function can be pr ovided in an explici t form (JNikolic et al 



2009) 



20111 ). Here we adopt the later 



or in an equivalent ordinary differential equation form as in (jTalathi et al. 
mathematical description as it provides a clear advantage in implementing various optostimulations protocols. 
In equation^ S o (0) = O.5(l + tanh(12O(0-O.l))) is a sigmoidal function and 9(t) = £\ ®(t-ti on )®(t-t ioff ) 
describes the optostimulation protocol; Q(x) = 1 if x > else 9(a?) = is the Heaviside function while ti on 
and ti of . are the onset, respectively the offset times of the ith optostimulation pulse. The constant tchr2 is 
the activation time of the ChR2 ion channel with typical values on the order of a few milliseconds. 



The ChR2 photocurrent is given by: 

-fchR.2 = V(giOi+g 2 o 2 ) 

= Vg\(oi +702); where 7 = — (10) 

The only parameter of the 4-state model described above that can be identified directly from the 
experimental data is G r = l/r r . Since the complexity of the 4-state model as represented in equa- 
tion [5] does not allow for estimation of closed form analytical solutions for the 4-state model parameters, 
a = {Pi, P2, Gdi, Gd2-, ei2, e2i, 7, gi, TchR.2}, we employ the following empirical estimation procedure to esti- 
mate these parameters. 

First, we establish an empirical formula to describe the photocurrent profile obtained from the experimen- 
tal data for the ChR2 photocurrent measurements obtained for a long term (Is) continuous optostimulation 
(I omp i) pulse and for a short term (2ms) optostimulation (I C mp2) pulse as follows: 

-ton)' 



/ t*-*on) \ 

4mpl(*) = ^pcak ( 1 - e -"so J 



®(t-t on )@((ton+t p )-t) 

(t-lto n +t p )) 



I pea . k \R+(l-R)e 

Q(t-(t on + t p ))@(t oS -t) 

Rlpe^e'^^Oit-tos) (11) 



^emp2\*V ^peak 



C*-*on? 

(1-e ^ ) (12) 



(t-Cw+tg)) 

0(i-i on )e(t p -i) + e ^ff e(i-tp) 



where is the Heaviside step function, i on is the time when the optostimulation pulse is turned on, t p is the 
time when ChR2 photocurrent reaches its peak value after the light is turned on and t g is the time when 
the optostimulation pulse is turned off. 

The first term in equation[TT]describes the rise of ChR2 photocurrent from zero to peak after the initiation 
of the optostimulation pulse. The second term simulates the decay of the ChR2 photocurrent from the peak 
to the steady state I p i a t and the last term accounts for the decay of the photocurrent from / p i a t to zero 



after the optostimulation pulse is turned off. Similarly, in equation 1121 the first term accounts for the rise 
of the photocurrent from zero and approaching peak value following brief optostimulation and the second 
term describes the exponential decay to zero from the maximum photocurrent value reached during the 
light on phase of optostimulation. Wc note that in writing equations [TT] and 1121 we have assumed that 
2 < Trise << 1000 ms, a valid assumption for all the ChR2 variants that we have studied here. 

We begin with an initial set of parameter values {a} (that obey the constraints detailed below) and 
numerically solve equation [§] to obtain IchR2(t,a) using equation 1101 with V = Vhoid for both long term and 
short term optostimulation protocols. We then evaluate the following objective function: 

C(a) = E 1 {a) + E 2 (a) + E 3 (a) (13) 

where: 



Ek(a) = 100* Ji- f l [Ic*R2(t,a) - I emp i(t)} 2 dt 
E 2 (a) = 100* Ji- / " [IchR 2 (t,a) -/cmp 2 (t)] 2 ^ 



(14) 



are the Root Means Square Deviations (RMSD) % errors of the photocurrent elicited by T\ — 1 s and T 2 = 2 
ms optostimulation pulse relative to their empirical profile and 



|-tpeak chR2 (QQ ^peak emp 
In 



E 3 (a) = 100 * * L1 "y cmpl (15) 

ipeak ompl 



is the percent error of the photocurrent peak achieved with prolonged optostimulation. 

We solve the constrained optimization problem: argmin C(a), using a global search method (genetic 

a 

algorithm) under the following constraints: 1) all parameters must be positive numbers, 2) for given light 
intensity (see table [!} , the values for the parameters P± and P 2 must not exceed a maximum value P ma x 
(Pi < P ma x and P 2 < P m ax) (evaluated using formula (38) from Appendix for maximum quantum efficiency, 
e = 1, and minimum degree of light absorbtion and scattering, wi oss = 1). 

The parameters are further optimized using a constrained local gradient descent algorithm (fmincon in 
Matlab). In Tabled we present the values for {a} obtained using the procedure described above for each 
of the ChR2 variants. 
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Experimental 
Data 


ChR2 
Variant 


Pi 

(ms _1 ) 


P 2 

(ms^ 1 ) 


Gdi 

(ms" 1 ) 


Gd2 

(ms- 1 ) 


ei2 

(ms -1 ) 


e2i 

(ms -1 ) 


(ms -1 ) 


TChR2 

(ms -1 ) 


7 


Si(*) 


Gunaydin 
(Gunavdin et al. 2010) 


ChRwt 
ChETA 


0.0641 
0.0661 


0.06102 
0.0641 


0.4558 
0.0102 


0.0704 
0.1510 


0.2044 
10.5128 


0.0090 
0.0050 


9.3458e-005 
le-003 


6.3152 
1.5855 


0.0305 
0.0141 


0.1136 
0.8759 


Berndt 
(Berndt et al. 2011') 


ChRwt 
ChET/TC 


0.1243 
0.1252 


0.0125 
0.0176 


0.0105 
0.0104 


0.1181 
0.1271 


4.3765 

16.1087 


1.6046 
1.0900 


9.3458c-005 
3.8462e-004 


0.504 
0.3615 


0.0157 
0.0179 


0.098 
0.5599 



Table 3: Parameters of the 4-state model. The parameters of the 4-state model have been determined for ChRwt and the 
two fast variants ChETA and ChRET/TC using experimental data (sec Table [TJ and search algorithms to best fit the profile 
photocurrent provided by equations II II and 1121 (*) This parameter is considered a variable in Section 3.2. 



2-4- Neuron Models 

The neural response to optostimulation is evaluated usin g two neuron models: 1) A single compartment 



fast spiking interneuron model of Wang and Buzsaki (WB) (jWang and Buz ski 



1996 ) and 2) A single com- 



Golomb et al 



2006 ). Each neuron 



partment hippocampal pyramidal neuron model of Golomb et al., (Gol) f( 

model with the addition of light dependent ChR2 ion channel current ichR2 can be described using a system 

of differential equations of the form: 



cV = Inc-I B (V,n)-IctiBa(V,B) 

ri = G(V, n) 

s = H(V,s) 

moo = U(V) 



(16) 



where V is the membrane voltage, n S [0, 1] and m^ G [0, 1] are the vectors of gating variables for the ion 
channels present on the neural membrane with hnite rise-time constant and instantaneous rise-time kinetics 
respectively; , s e [0, 1] is the vector of transition states of ChR2, I g is the sum of the membrane ion channel 
currents, IchR2 is the light activated ChR2 ion channel current, c is the membrane capacitance and Ir>c is 
the constant DC bias that controls the excitability of the neuron. 
For the WB model the membrane ion currents are, I g = {/Na,^K,-^L}, where 



Inu = gNam^hiV - E Na ) 
Ik = g K n\V-E k ) 
II = g L (V-E L ) 



(17) 
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4>(a h (l-h)-p h h) 



The gate variables are {n, h, m^}, where 

dh 

lit 

dn 

— = (p(a n {l - n) - p n n) 

Woo = f 1 -- (18) 

Otjn i Pm 

The functional form for the membrane voltage dependent rate constants {a m , (3 m ,ah, (3h,&n, ftn} are gives 
as: 



a m (V) 


-0.1(1/ + 35) 


(exp(-0.1(V + 35))-l) 


MV) 


= 4exp(-(V + 60)/18) 


a h (V) 


= 0.07(cxp(-(F + 58))/20) 


0h(V) 


= l/exp(-0.1(V + 28) + l) 


a n (V) 


-0.01(V + 34) 


(exp(-0.1(V + 34))-l) 


Pn(V) 


= 0.125cxp(-(F + 44)/80) 



(19) 

The model parameters are: Ej^ a = 55 mV, Ek = —90 mV, El = —65 mV, gNa = 35 mS/cm 2 , gx = 9 
mS/cm 2 , gi =0.1 mS/cm 2 , c = 1 /iF/cm 2 , = 5. For all the simulations, we set Il>c = — 0.51/iA/cm 2 , 
which ensures that the resting neuron membrane voltage is V Ics t = — 70 mV in absence of optostimulation. 
For the Gol model the membrane ion currents are, I g = {/Na, ^Kdr, Ih, I A, Im}, where 

iNa = gNam^hiV - E Na ) 

h = g L (V-E L ) 
iKdr = gpcdrn (V - Ek) 

iNaP = gNaPP 3 oc h ( V ~ E Na) 

I A = g A alb(V - E K ) 

hi = g M z(y-E K ) (20) 
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The gate variables are {n, h, b, z, moo,Pooi Qoo}, where 

h = <i> * (r(y,e h ,<T h ) -h)/{i.o + 7.5 *T{v,T th , -6.0)) 

b = {T(V,e b ,a b )-b)/T b 

n = <f>*(T(V,e n ,a n )-n)/(1.0 + 5.0*r(V,T tn ,-15.0)) 

z = (T{V,9 z ,a z )-z)/T z 

moo = T(V,6 m ,a m ) 

Poo = r(V, 6p,<7 p ) 

aoo = T{V,9 a ,a a ) (21) 

with T(V, 6, a) — A (v _ fl) x and the parameters: 9 m = —30; <7 m = 9.5; 9 p = —47; a v = 3; 9^ = 

1 + ex P(, ~5 j 

-45; ovi = -7; 5„ = -35; a n = 10; a = -50; a a = 20; fc = -80; afe = -6; 2 = -39; a z = 
5; n = 15; r 2 ; r t/l = -40.5; r t „ = -27. 

The remaining model parameters are: Ejss a = 55 mV; Ek = —90 mV; El = —70 mV; (f> = 
10; 9Na = 35 mS/cm 2 ; gjvaP = 0; g Mr = 6 mS/cm 2 ; g L = 0.05 mS/cm 2 ; g^ = 1.4mS/cm 2 ; g M = 
lmS/cm . For all the simulations with the Gol model, we set Inc = 0.12/iA/cm 2 . which ensures that the 
resting neuron membrane voltage is K-ost = — 70 mV in absence of optostimulation. 

Numerical Methods: All the simulations are performed using the forth order Rungc-Kutta method 
with a time step of 0.05 ms implemented in Matlab R2011b. The matlab codes will be made available for 
download at ModclDB. 

3. Results 

3.1. Modeling the response of non- excitable cells to optostimulation 

We begin by examining the extent to which the 3-state and the 4-state transition rate models can capture 
the main features of the ChR2 photocurrent elicited in a non-excitable (oocyte) cell by a prolonged (Is) 
optostimulation protocol. For the 3-state model, using the parameters obtained in Table [5] and assuming 
that all ChR2 molecules are in the dark adapted closed state C prior to the optostimulation, we find 
that, independent of the ChR2 variant, the 3-state model is unable to reproduce an important feature 
of experimentally measured ChR2 photocurrent, namely the plateau to peak ChR2 photocurrent ratio R 
(Figure (5J panels A2 and Figure 10 panel A2 from Appendix). Similar findings have been reported for 
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ChRwt photocurrent kine tics resulting from a 3-state model (JNikolic et al.l 120061) . In order to remedy this 



discrepancy, Nikolic et al., (JNikolic et al. 



2006) hypothesized that the recovery rate G r for the ChR2 molecules 



to transition from the desensitized state D to the dark adapted closed stat e C is a function of light intensity 



While this hypothesis is not implausible as stated in (INikolic et al 



20091 ) . to the best of our knowledge, it 



has not been yet supported by any experimental evidence. 
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Figure 2: Comparison example of the ChR2 photocurrent elicited by Is and 2 ms optostimulation. Al and 

Bl. Empirical data profile constructed for ChRwt(black) and ChETA(red) variant using equations ?? and ?? as well as the 
experimental data provided in Tablc[T]is displayed for a continuous Is optostimulation (Al) and for a brief 2 ms optostimulation 
(Bl) . The two normalized experimental photocurrent profiles are matching the experimental results reported in {Gunavdin et all 
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l2010h Figure 3c and Figure 2d. A2 and B2. Photocurrent generated by the 3-state model starting from ideal initial conditions 
((I): C(0) = 1; O(0) = 0; D(0) = 0. ) for ChRwt and ChETA variant for Is (A2) respectively 2ms (B2) optostimulation. A3 and 
B3. Photocurrent generated for ChRwt in comparison with ChETA, by the same 3-state model but starting from special initial 
conditions ((II): ChRwt: O(0) = 0.0132; O(0) = 0.0023; D(0) = 0.9845; ChETA: C*(0) = 0.0251; O(0) = 0.0085; D(0) = 0.9664.) 
evaluated in Appendix, Section 6.3. A4 and B4 Photocurrent generated for the same ChR.2 variants using the 4-state model 
with ideal initial conditions ((I): Oi(0) = l;Oi(0) = 0; O2 = 0.;C2(0) = 0. ) ). For a similar example comparison of ChR2 
photocurrent elicited by ChRwt and ChRET/TC see Figure 10 in Appendix. 



We next investigate whether the feature R, of ChR2 photocurrent, can be captured using a different set 
of initial conditions. We note that, the dependence of the 3-state model solution on the ini tial condition 



of th e state variables was not evaluated in the original formulation of the 3-state model by (JNikolic et al 



2006). In order to address this question, we analytically evaluate the set of initial conditions required to 
correctly reproduce the ratio R (see Section 6.3 of Appendix). We find that the initial conditions that 
generate experimentally observed ChR2 photocurrents with reproducible features, correspond to the case 
when all the ChR2 molecules are in the desensitized state D at the beginning of the optostimulation protocol, 
suggesting the possibility that the experiments might have been carried out under non-ideal conditions, i.e, 
all ChR2 molecules have not relaxed back to the closed state before the prolonged optostimulation protocol 
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was run (Figure [5J panels A3 and Figure 10, panel A3 from Appendix). 

In order to determine whether this is indeed the case, we conduct the following simulation experiment: 
we begin with initial conditions for the 3-state model in the desensitized state D. We then present two 
prolonged optostimulation pulses of 1 s in duration each, separated by time < At < 50 s. We then 
calculate the % recovery Rec(At) = 100 * AI2/AI1 where AZi^ are the differences between the peak and 
steady state photocurrents for the two optostimulation pulses respectively. In Figure |3] A and B we present 
these results in comparison with the corresponding results obtained for simulations started with ideal initial 
conditions (see Figure [3] C and D). It is clear from Figure [3J3 that the photocurrent profile produced by 
the 3-state model beginning with all ChR2 molecules in the desensitized state D exhibits a recovery percent 
grater than 100 (Rec(At) > 100 for At > 0) . When the 3-state model is simulated with the special initial 
condition of all ChR2 molecules being in the desensitized state, at the end of first optostimulation protocol, 
more ChR2 molecules are present in the closed state and ready to maximally respond to an impending 
optostimulation protocol res ulting in Rec(At) > 10 0. This finding is in stark contrast with the experimental 



results of Gunaydin et al., (jGunavdin et al 



2010I ). The experimentally observed recovery curve for ChR2 
photocurrents can be reproduced if we assume that all the ChR2 molecules at the beginning of the experiment 
arc in the dark adapted closed state C as shown in Figure [3P and[3p. 

Following from these observations, we conclude that while the 3-state model can exhibit ChR2 pho- 
tocurrent characteristics for a single prolonged optostimulation pulse under special circumstances, it is not a 
suitable model to mimic the transient phot ocurrent dynamics in response to multiple optostimulation pulses. 



We note that the key reason presented in (jNikolic et al 



20091) for the unsuitability of 3-state model was low 



model complexity with insufficient model parameters to capture all the essential features of ChR2 photocur- 
rent profile, specifically, the bi-exponential photocurrent decay. Our findings on the other hand suggest 
that for all the ChR2 variants analyzed here, in terms of model complexity the 3-statc model is sufficient 
to capture all the essential features of ChR2 photocurrent profile. However, this requires non-equilibrium 
initial distribution of ChR2 molecules across the 3-states in dark adapted conditions. Further support for 
our claim of the unsuitability of 3-state model will be presented in Section 3.2. 

In the case of the 4 state model, we arc able to employ the strategy described in Section 2.3 to identify 
the set of model parameters that yield ChR2 photocurrents matching the experimental recordings for an 
ideal set of initial conditions i.e., all the ChR2 molecules are in the dark adapted closed state C\ satisfying 
the initial conditions c\ = 1 and o\ = 0-2 = C2 = 0. The following general conclusions can be drawn from 
the estimated 4-state model parameters (see Figure 0] for a visual comparison of these values): 1) For all the 
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Figure 3: Example of recovery curves in the 3 state model A. and C. Recovery curves have been evaluated for ChRwt 
(black) and ChETA (red) for special initial conditions (A) respectively for ideal initial conditions (C) at the beginning of the 
optostimulation protocol. B Example of photocurrent curve obtained for ChRwt when the protein is in special initial condition 
at the beginning of the stimulation and two pulses of 1 second each are applied with a dark period of 1 s in between. The 
recovery (Rec(At) = 100 * AZ2/A/1) is inconsistent with the experimental observations. D. Similar example of photocurrent 
curve as in B obtained when the protein is in ideal initial conditions at the beginning of stimulation. The recovery is consistent 
with the experimental observations but the ratio R of photocurrent peak to steady state is not satisfied. 

variants P\ > P2, which underscores the idea that ChR2 molecules exhibit two distinct dark-light adapted 
closed states with the closed state C\ being more light sensitive than the closed state C2 to optostimulation 



as was suggested by Nikolic et al., (JNikolic et al 



2003 ). 2) The fast transition of ChR2 photocurrents from 



the pe ak value to steady stat e for both variants relative to ChRwt o bserved in the experiments of Gunaydin 



et al., (jGunavdin et al. 



20101 ) and Berndt et al. (JBerndt et al 



20111 ) are explained via significant differences 



in the transition rate parameters Gdi and/or ei2 across the two variants. Specifically, Gj™ wt >> G^ ETA , 
e ChRwt << e ChETA an( j e ChRwt << e ^ BET ' TC ^ indicating that the transition from the dark adapted open 
state 0\ (which accounts for the peak ChR2 photocurrent) into the light adapted open state O2 (which 
accounts for the steady state ChR2 photocurrent) is the preferred mode of operation for the ChETA variant 
over the situation when the dark adapted open state transitions to the dark adapted closed state. 3) Both 
variants exhibit higher conductance (quantified here by the parameter g{) relative to ChRwt. While in 
the case of ChETA the expected increase in the maximum amplitude of the photocurrent is suppressed 
by an enhanced value of the parameter Gdi (controlling the decay of the second open state), in the case 
of the ChRET/TC variant this suppression does not occur (as 67 ?j lRwt ~ G W9 ) which explains the 



significantly enhanced photocurrent peak reported for this variant ( 



Berndt et al. 



20111 ). Finally, 4) we notice 



that both variants show decreased values of the time constant of the channel opening relative to ChRwt 



16 



( T ShR2 A ^ T cftfl2°' anc ^ T chR2 < T ShR2 t ) suggesting that the time scale of the conformational changes 

leading to the channel opening is significantly faster in the newly developed variants relative to ChRwt. 
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Figure 4: Comparison between the coefficients of the 4 state model. A. Comparison between the 4-state model 
coefficients evaluated for ChETA (red) and ChRwt (black). B. Similar comparison between the coefficients of ChRwt (black) 
and the fast variant ChRET/TC (blue). Because of significant differences among the values of the parameters, two scales have 
been used to visualize their comparison. In both panels, the parameters to the right of the dashed line obeys the scale on 
the right hand side while the parameters to the left of the dashed line obeys the left hand side scale; difference between the 
corresponding parameter values for ChRwt and the variant (ChETA or ChRET/TC) larger than 50% but smaller than 3 times 
the value are indicated with one star and differences larger than 3 times the value are indicated by two stars. Actual values for 
all parameters are reported in Table [3] 

It is important to note that for the ChR2 photocurrent profile generated using either 3-state and 4-state 
models, the parameter g±, expressing the overall conductance of the ChR2 variant, controls the maximum 
value of the photocurrent amplitude and is independent of other characteristics of the photocurrent dynamics 
governed by the transition rate models such as the activation (r r jse)j inactivation (tj;„) and deactivation (t d //) 
time constants and the ratio R of peak to steady photocurrent amplitude. This parameter therefore, can 
capture the degree of expression of the ChR2 variant in the membrane of the neural cell and is assumed to 
be variable across experiments. 

3.2. Modeling the response of excitable cells to opto stimulation 

We now investigate whether the two transition rat e models can capture the experimental re sults on the 



neural response to various optostimulation protocols (jGunavdin et al 



2010; 



Berndt et al 



20111) . We begin 



our inv estigation by incorporat ing the 3- and 4-state ChR2 photocurrent transition rate models into the WB 



model (jWang and Buzski 



1996}). In Figure [SJ we compare our simulation results with experimental data on 



neural activity elicited in the fast spiking hippocampal inter neuron models express ing ChRwt and ChETA 



for various short duration (<<1 s) optostimulation protocols (jGunavdin et al 



20101 ) . We make the following 



observations: for low frequency optostimulation protocols (< 20 Hz) we find that both the 3- and 4-state 
models are able generate neural response from the WB model neuron consistent with the experimental 
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recordings. For instance, for 10 Hz optostimulation frequency, the simulated neural response using the 4- 
state model for ChR2 kinetics (see Fig. [SJA. and [5j3 the top pan els) matches very well the experimental 



20101 )'). Similar results 



observations for the same optostimulation protocol (Figure 4a in (jGunavdin et al. 

are obtained when the 3-state model is employed to model the ChR2 photocurrent kinetics (see top panels 

in Figure [5p and[5p). 

For higher frequency (> 40 Hz) optostimulation protocols, significant differences in the neural activity 
elicited by the two transition rate models are revealed. In particular, the 3-state model is not able to capture 
the bursting features observed in the neural activity of interneurons expressing ChRwt in response to 80 
Hz optostimulation protocol (se e Figure \5Jp middle pa nel in comparison with Figure [5j\ middle panel and 



Figure 4a middle-left panel in (jGunavdin et al. 



20101 )); neither is the 3-state model able to capture the 



neural activity of interneurons expressing ChETA. For the later case, the neural response exhibits features 
inconsistent with the experimental findings including extra spikes at the beginning of the stimulation protocol 
and missed spikes at the end of the stimulation protocol (see Fig. [5p middle panel). These inconsistencies 
are not observed when the 4-state model is i mplemented in the WB model neuron (Figure [3)3 middle panel 



compared to Figure 4a middle-right panel in (|Gunavdin et al 



20101) ). 



For very high frequency (> 150 Hz) optostimulation protocols, the 4-state model is again able to better 
capture the features of neural response for both the ChR wt and ChETA mutan t (Figure EJA and [5)3 bottom 



2010I )). For instance, the plateau 



panels in comparison with Figure 4a bottom panel from (IGunavdin et al. 
potential observed for neurons expressing ChRwt is well reproduced by our simulation results produced 
using the 4-state model (see Figure [5j\ bottom panel) and at the same time our simulation experiments 
match the experimental findings from the neuron expressing ChETA in the sense that the plateau behavior 
is absent (Figure [5)3 bottom panel) . On the other hand, the simulated neural response of the WB model 
neuron with a 3-state model implementation of the ChRwt and ChETA variant shows several artifacts such 
as, a decrease of the plateau potential when a larger number of stimulation pulses is presented to the neuron 
expressing ChRwt, depolarization block at the beginning of the stimulation protocol and inconsistencies in 
the interspike intervals observed in neurons expressing ChETA (Figure [5p and [5)3 bottom panels). These 
inconsistencies are primarily driven by the limitations of the 3-state model to mimic the recovery rate in 
response to multiple short duration optostimulation light pulses, as explained in the Section 3.1. 

We have also investigated the neural response of the WB model neuron to prolonged (more than 50 
stimuli) optostimulation protocols. We observe that both the 3-state and the 4-state model can capture the 
behavior of spike failure observed in experimental data on interneurons expressing ChRwt (Figure [5] panels 
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A and C in comparison with the experimental results presented in (jGunavdin et al.ll2010l ). Figure 3e) and 
the correct neural response observed when the fast variant ChETA is expressed (see Figure [BJ3 and|Bp). 
However, with the 3-state model implementation, we observed a transient spike adaptation phenomenon 
in the spiking response of the WB model neuron at the beginning of the optostimulation protocol (more 
pronounced for 3-state model implementation of ChRwt), a phenomenon that has not been reported in 



experimental recordings (IGunavdin et al. 



2010) 



We have also conducted a comparative study to evaluate the performance of the 3-state and 4-state mod- 
els in mimicking experimental data on neural response of hippocampal pyra midal cells expressin g ChRwt 



and a fast variant labeled ChRET/TC to various optostimulation protocols (JBerndt et al 



2011). We in- 



corporated a 3-state and 4-stat e model for ChRwt and the variant ChRET/TC into the Col model for 



hippocampal pyramidal neurons (|Golomb et al. 



J2006J ). We find that when the 3-state model for both ChRwt 
and ChRET/TC is implemented in the Col model neuron, the model neuron elicits a significant number of 
additional spikes at the beginning of a 40 Hz optostimulation protocol, a feature that is inconsistent with 
experimental findings. On the other hand, the 4-state model implementation produces neural response that 
closel y matches those ob served experimentally (Figure [7R and 03 respectively in comparison to Figure 4A 



from ( Berndt et al 



l20iin ). 

It is evident from the number of examples presented above that the 4-state model is better able to account 
for the multitude of neural response features elicited by optostimulation protocols in both excitatory and 
inhibitory model neurons. We therefore focus on the 4-state model and further investigate the extent to which 
several other features of neural activity generated by optostimulation can be captured using the 4-state model. 
We evaluate whether the number of extra spikes elicited by optostimulation protocols of different frequencies 
when 4- state models for C hRwt and ChETA vari ant are expressed in WB model interneuron match those 



reported in experiments of (jGunavdin et al 



20101) . We find a good correspondence with t he experimental 



findin gs for both the ChRwt and ChETA variants (Figure [8]A. in comparison to Figure 4b in (jGunavdin et al 



2010h ). We also report on the ability of the 4-state models for ChRwt and ChETA to mimic experimental 
findings on the plateau potential elicited by the WB model neuron in response to optostimulation with 
various frequencies (F igure [BP). As reported in the experimental findings of Gunaydin et al., (Figure 4c in 



(JGunavdin et al. 



20101 )). we find that the WB model neuron expressing ChRwt consistently generates higher 



values for the plateau potential as compared to those generated by the WB model neuron expressing ChETA. 



Finally, we report results on the 4-state model's ability to capture the dependence of the firing success 
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Figure 5: Comparison between the simulated interneural response for different optostimulation protocols when 
4-state respectively 3-state model accounts for ChR2 dynamics. A. Simulated neural response to 3 (upper panel), 20 
(middle panel) and 40 (lower panel) stimuli of 2 ms each presented at 10, 80 and 200 Hz optostimulation frequency obtained 
when the 4-state model is employed to account for ChRwt dynamics; the overall conductance was g\ = 100 (upper panel), 
gi = 40 (middle panel), respectively g\ = 20 (lower panel). B. The response to the same optostimulation protocol is simulated 
in intcrncuron model expressing ChETA; the over all conductance was g\ = 70 for all simulations. These results are in excellent 
agreement with experimental data presented in jGunavdin et aU 120101) Figure 4 a). C. and D. Results obtained for the 
same optostimulation protocol when the 3-state model is representing the ChRwt (C) and ChETA (D) kinetics. The overall 
conductance has the following values: in C ji = 4.8 (upper and middle panel), g± = 4 (lower panel); in D ji = 1.2 (upper 
panel), gi = 2.2 (middle panel) and gi = 4 (lower panel). The differences between the results obtained for the 3 and 4- state 
model arc discussed in detail in the text. 
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Figure 6: Comparison of neural response to a prolonged optostimulation protocol. The neural activity elicited in 
interneurons expressing ChRwt (black) and ChETA (red) by a train of 60 optostimuli presented at 20Hz, each pulse of 2 ms 
duration, is presented in comparison when the 3 and 4 state models are employed to account for ChR2 kinetics. A and B 
Response of interneurons expressing the two types of channclrhodopsins when the 4-state model is employed to mimic the 
ChR2 kinetics. The overall con ductance was qi = 10 .8 (A) and gi = 70 (B). These results are in good agremcnt with the 
experimental data presented din l lGunavdin et aUl2010f) Figure 3 e). C and D Similar results obtained when the 3-state model 
is used to account for ChR2 dynamics. The overall conductance was g\ = 4 (C) respectively g\ = 2.2 (D). 
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Figure 7: Comparison of pyramidal response to a prolonged optostimulation protocol. The simulated neural response 
to an optostimulation protocol comprising 60 stimuli of 2 ms each presented at 40 Hz is evaluate for different configurations. A 
and B Results obtained in pyramidal neuronal model by employing the 4-state model to account for ChRwt (A) and ChRET/TC 
(B) kinetics. The o verall conductance w as g\ = 4.8 (A) and g\ = 33 (B). These results are consistent with the experimental 
findings reported in l|Berndt et al.ll201lf l. Figure 4 A. C and D Similar results obtained in pyramidal neurons when the 3-statc 
model has been employed to account for ChR2 dynamics. The overall conductance was g\ = 12.6 (C) and g\ = 12 (D). 
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Figure 8: Evaluation of the number of extra spikes and plateau potential. A. The number of extra spikes have been 
evaluated when optostimulation trains of 40 stimuli of 2 ms each and various frequencies have been applied on interneurons 
expressing ChRwt(black) and ChETA(red). In ChRwt simulation s presented here tckri = 0.5. The results are comparable 
with the experimental findings presented in llGunavdin et aUl2010h Figure 4b B. Modeling results when the plateau potential 
is eva luated in the same conditions as in A. are qualitatively similar with the experimental results presented in (jGunaydin ct al. 
l2010h Figure 4c. 



rate of a pyramidal neuron on the light intensity of optostimulation protocol (see Figure |H|) . We find that 
the general trend, observed experimentally, of decrease in the success rate of optostimulation protocol to 
induce an action potential spike in pyramidal neuron is captured well in our simulation using the 4-state 
model for ChRwt and ChRET/TC incorporated in the Gol model neuron. However, we also find certain 
discrepancies. For example, for higher optostimulation frequencies (>60 Hz), for all light intensities tested, 
the firing success rate of the Gol model neuron using e ither of the two Ch R2 variants is consistently higher 

20111 ) ). While, for low intensity (6.7 



than those reported in experiments (Figure 4B from (jBerndt et al. 
mW/mm 2 ) and low frequency (< 30 Hz) optostimulation protocol, the firing success rate observed in Gol 
model neuron is lower than those reported in experiments. These inconsistencies may be due to processes 
such as the dependence of quantum efficiency of the ChR2 protein on light intensity (e = e(/)) or variability 
in the degree of absorption and diffraction as a function of light intensity intensity (wi oss = wi oss (I)), which 
are not captured by the 4-state model employed here. 



3.3. Degeneracy in the 4-state model 

Our analysis thus far suggests that the 4-state model for ChR2 photocurrent kinetics is better suited to 
account for the dynamical features of the experimentally observed neural response to various optostimulation 
protocols. This model gives rise to a general solution with three time constants in the light on condition 
and two time constants in the light off condition. We note however that the majority of experimental 
studies concerned with the molecular kinetics of light activated ion channels report a single time constant 
for the decay of the ChR2 photocurrent from peak to steady state when light is turned on and a single time 
constant for the decay of the photocurrent from steady state to zero when light is turned off. These data 
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Figure 9: Modeling results of neural response to different optostimulation frequencies with light of different 
intensities. Modeling results of firing success rate simulated in pyramidal neurons expressing ChRwt (black) and ChRET/TC 
(blue) using the 4-state model and optostimulation protocols identical of 60 light puls es of 2 ms each and different stimulations 
frequencies. The results are comparable with the experimental findings presented in l|Berndt et al.ll201ll) Fig. 4B. 

suggest that a 3-state model is an appropriate mathematical construct to mimic ChR2 photocurrcnt kinetics. 
However, as we have reported above, the 3-state model is unable to capture all the feature of neural response 
to optostimulation protocols. We also note that a most generic 3-state model than considered here (not 
reported here) also fails to capture many of the key characteristics of neural response to optostimulation. 
The 4-state model can exhibit mono-exponential decay characteristics for the ChR2 photocurrent kinetics 
under the following three special conditions: (1) The solution to the 4-state model exhibits degenerate time 
constants and/or (2) the amplitude of the currents associated with one of the time constants is negligible 
and/or (3) there is at least an order of magnitude difference in the time constants associated with the solution 
of the 4-state model corresponding to a rapidly decaying current component, that is not captured in the 
empirical fit to the experimental data. 

In order to evaluate whether any of these conditions are satisfied for the set of parameters reported 
in Table [3] for the 4-state model, we obtain an approximate analytical solution to equation [9] under the 
assumption that the conformational changes mediating the channel opening are considered instantaneous, 
i.e., s = 1. This approximation holds as long as the time constant tckri is not significantly larger than the 
peak time of the photocurrent t p , a condition satisfied by three of the ChR2 variants invest igated in this 



paper (ChRwt and ChRET/TC model ed based on the data 



arovided in 



modeled based on the data provided in 



Gunavdin et al 



Berndt et al. 



|201ll) and ChETA 



(|2010l )). We begin by identifying a semi- analytical 



solution for the light on condition in the 4-state model (Appendix, Section 6.5 for details). The photocurrent 
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can be expressed as: 

Ion(t) = / r isc(o„)e- A — <°"' f + /s(o„)e- A (-'' 

+ /f(on)e- Af< °">*+/plat (22) 

where, T r w on ) = AT , . represents the time constant associated with the rising phase of the photocurrent 
(also labeled as the activation time constant), r s ( on ) = A~, n \,Tf( on ) = A« ., are the slow and fast time 
constants associate with the bi-exponential photocurrent decay from peak to plateau (also referred to as the 
inactivation time constants), / r ise(on) is the amplitude of the rising phase of the photocurrent, I s ( on ) is the 
amplitude of the slow component of the photocurrent decay from peak to plateau, /f( on ) is the amplitude 
of the fast component of the photocurrent decay from peak to plateau and I p i a t is the steady state plateau 
photocurrent. The expression for all these photocurrent components can be found in Section 6.4 of the 
Appendix. 



2009 ) which 



For the light off condition we use the solution provided in Nikolic et al. 2009 (jNikolic et al 
reads: 



/off - I s{ oS)e- K ^ t + / f (off)e- Af <° ff)t ; (23) 

where T a t g) = A", « , Tft g) = A~ off , are the slow and fast eigenvalues associate with the bi-exponential 
photocurrent decay from plateau to zero, / s ( ff) and if( ff) are the amplitudes of the slow and fast photocurrent 
decay from plateau to zero after the light is turned off. For both solutions the analytical expressions of the 
time constants and the current amplitudes can be found in Appendix, Section 6.5. 

In Table 21 we present the values of the fast and slow time constants and the corresponding amplitude 
values for the photocurrent decay from peak to plateau under the light on condition and from plateau to 
zero under the light off condition. We find that for each variant, for light off condition, the time decay of the 
slow component (r s = A J 1 ) provided by the analytical solution is very similar in value with the one reported 
experimentally; the second (usually unreported) time decay (ry = AT ) is considerably smaller, implying 
that its associated photocurrent component decays on a much faster time scale relative to the experimental 
measurements. Furthermore, the amplitude of this photocurrent component is also significantly smaller as 
compared to that of slower component of the photocurrent (sec Table 0] for a comparison). 

For light on condition ChRwt and ChRET/TC variant exhibit degeneracy, namely the fast and slow time 
constants are of same order in magnitude and their sum can be well- approximated using a single exponential 
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Light On 


Light Off 


ChRwt 


T in = 9.6 (experimental) 


T ff = 11.1 (experimental) 




t s = 10.91; A s = 0.736 
T f = 7.47; A f = -0.756 


t s = 11.255; A s = 0.036 
T f =0.166; A f = -0.0004 


ChETA 


Tin = 15 (experimental) 


T ff = 5.2 (experimental) 




t s = 14.91; A s = 0.0065 
r/ = 4.65; ,4/ = -0.0046 


t s = 6.625; A s = 0.0043 
Tf =0.095; Af =0.0001 


ChRET/TC 


Ti n = 11 (experimental) 


T ff = 8.1 (experimental) 




t s =8. 11; ^ =0.5457 
r/ = 7.17; A f = -0.5494 


t s =8.36;A S = 0.0105 

Tf =0.06; Af = -3.3e-005 



Table 4: Comparison between the experimental and analytical time decay constants. The photocurrent decay from 
peak to plateau Ti n and from plateau to zero (t //) for each ChR2 variant is compared with the fast (Tf) and slow (ts time 
constants)derived analytically. The maximum amplitude of the photocurrent component associated with each time decay is 
also provided. 



decay function. For the variant ChETA no degeneracy is observed; the slow time constant is very similar 
with the experimentally observed time decay while the fast photocurrent component has a significantly lower 
time constant and amplitude (conditions similar with the ones observed for light off condition). 

Together, these results provide clues to the absence of reports on the bi-exponcntial decay time constants 
for the ChR2 photocurrent profiles in majority of in vitro experiments and at the same time provide support 
for a 4-state model as an appropriate model to mimic ChR2 photocurrent kinetics of various types of light 
activated channels. 

4. Discussion 



In this paper we have systematically analyzed the 3-state and 4-state transition rate models for several 
ChR2 variants. Using experimental data available from recently published experimental studies we first 
identified the 3-state and 4-state model parameters. We incorporate these models into model neurons and 
investigate the extent to which the dynamical features of model neuron response match to those observed 
experimentally for various optostimulation protocols. We demonstrate that the 3- state model present s sig- 



20061 ). but 



nificant limitations, not because of low model complexity as was reported earlier (jNikolic et al. 
rather due to its inability to mimic the correct ratio R of steady state to peak photocurrent in the presence of 
prolonged light stimulation and as a result its failure to account for for certain dynamical features of neural 
activity elicited by optostimulation of hippocampal excitatory and inhibitory neurons. In comparison, the 
implementation of the 4-state model produced better match to the experimental data. 

Two different strategics were adopted to estimate the parameters for 3-state and 4-state models from the 
available experimental data of ChR2 photocurrent profiles. The low complexity of 3-state model afforded 
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analytical expressions to uniquely determine the 3-state model parameters in terms of experimental data. 
On the other hand, the increased complexity of 4-state model warranted empirical procedure to estimate 
all the model parameters. Motivated by the overall better performance of the 4-state model, we performed 
further investigations to determine the conditions under which a characteristic specific to the 3-state model, 
namely the mono-exponential decay of the ChR2 photocurrent kinetics (often reported in the experimental 
literature) can occur in the 4-state model. 

For the variants investigated in this paper, there are two possible ways in which the 4-state model can 
account for the observed mono-exponential decay characteristics of ChR2 photocurrents: first, a degeneracy 
in the light on condition for two of the variants (ChRwt and ChRET/TC); second, or alternatively, a 
significant difference in the two decay time constants with one decay time constant matching those reported 
experimentally while the second time constant as well as the associated amplitude component being small 
enough to remain unidentified in the empirical fit to the experimental photocurrent profile. Together, these 
results suggest that the 4-state model is better suited to capture the complexity of the photocurrent kinetics of 
all ChR2 variants tested, possibly representing a generic mathematical framework to model the photocurrent 
dynamics of ChR2 proteins. 

Our detailed computational analysis suggests that some neural dynamical features elicited by optostimu- 
lation depend on the degree of expression of the ChR2 variant in the cell. This dependence is more prominent 
for ChR2wt and can be observed for optostimulation protocols performed on interneurons at rest (see Figure 
[5]). These observations suggest that in addition to the parameters accounting for the ChR2 kinetics (such as 
the activation and inactivation time constants, etc) the degree of ChR2 expression should be also considered 
a parameter of significant importance in computational models of neural activity elicited by optostimulation. 
It is also important to note that while employing the 4-state model provides results in excellent agreement 
with the experimental recordings across a multitude of optostimulation protocols, not all the features can 
be captured with a unique set of parameters. This aspect is particularly important for future optogenetic 
control strategies of neural function which generally rely on the consistency between the experimentally 
measured neural activity and the model based, predicted neural response. 

Progress in the field of genetics and molecular biology has resulted in the identification of a myriad of 
light sensitive proteins with enhanced kinetic features and improved sensitivity to certain light wavelengths. 
For the purpose of this paper, we have chosen to discuss two fast variants (ChETA and ChRET/TC) in 
comparison with the more common ChRwt for the following reasons: 1) the experimental studies describing 
the response of these proteins to optostimulation provided all the necessary data to uniquely determine 
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the parameters for the 3-state model and constrain the parametric search in the 4-state model; 2) these 
studies have also provided detailed information on the neural response (of two important classes of neurons, 
namely hippocampal pyramidal cells and PV-interneurons) to a multitude of optostimulation protocols, 
providing strong experimental datasets necessary for evaluating the efficiency of different model for ChR2 
photocurrent kinetics. Evidently, the modeling methods outlined in this paper and the analysis thereof, can 
be easily implemented in the case of any other variant of potential interest in research applications. 

Our results demonstrate the ability of 4-state transition rate model to mimic experimental findings 
across a multitude of optostimulation protocols applied on interneurons as well as pyramidal cells expressing 
different types of ChR2 variants. The ability to successfully bridge the experimental and computational 
observations is critical for future development of optostimulation control protocols. In this context, this 
paper, speaks in favor of the exciting possibility of using computational models to support and enhance 
the design of novel control strategies and optostimulation protocols to regulate normal and abnormal neural 
network activity. 
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6. Appendix 

6.1. Analytical solution of the 3-state model. 
The equations describing the model are: 



6 = P(l-o-d)-G d o (24) 

d = GdO — G r d 



The photocurrent is: 



I = V 9l o (25) 
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The equivalent theoretical solution is: 

I = Cie- Alt + C 2 e- X2t + I plat (26) 

where 

Ai = a - P; \ 2 = a + f3 (27) 



P + G d + G r (P-G d -G r )* P(X 1+ X 2 -P-G d ) 
Oi= ; P = y o G d L, r \ ipiat = t-t (28) 



and 



r DoP^i + (Ag - P - G d )(P - OqAQ 

Cl = A^-A,) (29) 

r n r P(P + G d - Ai-A 2 ) 



A1A2 

When ideal initial conditions arc satisfied (Dq = 0; Oo = 0) the solution presented in 
is recovered. 

6.2. Derivation of formula |3|) from the main manuscript. 

Wc start with the expression of the first eigenvalue given by equation ?? above: 



Nikolic et al 
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. P + G d + G r y/{P -G d - G r f - AG d G r 

Ai = ~ 2 ~ 2 (30) 



We multiply the equation above with 2 and rearrange the terms to obtain: 



P + Gd + G r - 2Ai = y/(P -G d - G r ) 2 - AG d G r (31) 

we raise the equation above to the second power and rearrange the terms to obtain: 

A? - Xi(G d + G r ) + P(G d + G r - Ai) + G d G r = 0; (32) 

which finally gives: 
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^ = A 1+ . G f d (33) 

Aj — U r — Ud 

6.3. Evaluation of necessary initial conditions for the 3- state model. 

We can find the initial conditions necessary to obtain a photocurrcnt which will exhibit the appropriate 
Ipeak I Ipiat ratio by use of the following conditions: 

1. We first find the coefficients of the homogeneous solution (Ci, C 2 ) which will satisfy the experimental 
data by solving the following system of equations: 



AiCie- Al *" + X 2 C 2 e- X2tp = (34) 

C e -Ait p q e -A 2 t p ]_ 

— h — hi = — 

Iplat Iplat R 

where the first equation represents the condition that the derivative of the photocurrcnt function is zero 
for the / = Ip ea k and the second equation instantiate the condition that the ration Ipeak / Iplat must match 
the ratio -5 provided by the experimental data. 

2. With C\ and C 2 determined above we can find the initial conditions (00,00) by solving the equations 
7 and 8 given in the previous section; then Cq — 1 — oq — do . 

The solutions for the all of the above equations have been evaluated symbolically in Matlab by using 
the function solve and then numerically by allowing the parameters to take the appropriate values in the 
symbolic solution. 

6.4. Experimental Results - Additional Information. 
6.4.I. Evaluation of the activation time constant T r i se . 

The evaluation of the time constant of the rising phase of the photocurrent from zero to peak has been 
performed by approximating the photocurrcnt curve with a mono-exponential function. Thus, wc can write: 

I(t)=I ma xO--e'^) (35) 

when the the photocurrent reaches maximum (at t = t p ) we can approximate: 
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0.99999; 



(36) 



which leads to: 



Trise — 



Zn(O.OOOOl) 



(37) 



6.4-2. Comparison between the photocurrent induced by continuous Is and brief 2 ms opto stimulation in cell 
expressing ChRwt and the fast ChRET/TC variant. 
We present in FigJT0]a comparison between the ChR wt and ChRET/TC photocurrent elicited by Is and 
2 ms continuous optostimulation. 
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Figure 10: Comparison example of the ChR2 photocurrent elicited by Is and 2 ms optostimulation. Al and 

Bl. Empirical data profile constructed for ChRwt and ChRET/TC variant using equations 11 and 12 from the main paper as 
well as the experimental data provided in Tablel is displayed for a continuous Is optostimulation (Al) and for a brief 2 ms 
op tostimulation (Bl). Th e two normalized experimental photocurrent profiles are matching the experimental results reported 
in lGunavdin et al.l J2010T> Fig 3c and Fig 2d. A2 and B2. Photocurrent generated by the 3-state model starting from ideal 
initial conditions ((I): C(0) = l;O(0) = 0;D(0) = 0. ) for ChRwt and ChRET/TC variant for Is (A2) respectively 2ms 
(B2) optostimulation. A3 and B3. Photocurrent generated for ChRwt in comparison with ChR ET/TC, by the same 3-statc 
model but starting from special initial conditions ((II): ChRwt: C*(0) = 0.0041; O(0) = 0.0037; D(0) = 0.9922; ChRET/TC: 
C(0) = 0.00156; O(0) = 0.0098; D(0) = 0.9746.) evaluated in Appendix, Section 6.3. A4 and B4 Photocurrent generated for 
the same ChR2 variants using the 4-state model with ideal initial conditions ((I): Ci(0) = l;Oi(0) = 0; O2 = 0.;C2(0) = 0. ) ). 



6.4-3. Dependence between the excitation rate (P) and light intensity (I). 



n r-, Vret(f> ea ret \ T 

P = er = e = —1 

Wl oss Wi oss hc 



(38) 
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where e is the qu a ntum efficiency of photon absorbtion (a typical value for rhodopsin is e ~ 0.5 



Braun and Hegemannl (|1999h ). F = <r ret (j)/wi oss is the number of photons absorbed by C hR2 molecule 
per unit time, a re t is the retinal cross-section (a re t ^ 1.2 x 10~ 20 m 2 



Hegemann et al 



(120051) ). wi oss is the 



measure of the loss of incidental photons due to scattering and absorbtion phenomena, <f> = XI /he is the 
photon flux per unit area, A ~ 480nm is the wave length of the light used in the stimulation protocol, 
I(mW/mm 2 ) is the light intensity, h — 6.626 x 10 -34 Js is the Planck's constant and c = 3x 10 8 m/s is the 
speed of light in vacuum. 

6.5. Derivation of the semi- analytical solution for the light on condition in the J^-state model. 
The 4 state model can be written as follows: 



o'i = Pi(l - c 2 - o x - o 2 ) - (G d i +ei 2 )oi + e 2 \o 2 
o' 2 = P 2 c 2 + ei 2 oi - (G d 2 + e 2 i)o 2 
C 2 = Gd202 - (P2 + G r )c 2 



(39) 



With the notation y = \o\ o 2 c 2 ] T , equation 1391 can be than expressed as follows: 



V = 



-(Pi+G dl + e 12 ) 621 -Pi -Pi 

ei2 -(G d2 + e 2 i) P 2 

G d2 -(P 2 + G r ) 



Pi 





= Ay + P 



The general solution that needs to be evaluated can be written as: 



y = Vc + y P 



(40) 



(41) 



where y c of the homogeneous (complementary) solution and y p is the particular solution of the system 
of equations (18). In the following we will evaluate both components. 

6.5.1. Finding the eigenvalues. 
The characteristic equation is: 
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det(A - XI) = 



(42) 



or 



which leads to 



■(Pi + Gdi + ei 2 )-A eai-Pi -Pi 

ei2 -(Gd2 + e 2 i) - A P 2 

G d2 -(P2 + G r )-\ 



= 



(43) 



- [(Pi + G dl + e 12 + \){G d2 + eai + A)] - Pie 12 G d2 + 
P 2 G d2 {P 1 + G dl + eia + A) + e 12 (e 21 - Px)^ + G r + A) = 



(44) 



This equation is solved symbolically in Matlab using the commend solve which gives the expressions 
for the solutions: Ai, A 2 , A3. The actual expressions are very elaborated, therefore they will not be included 
here. The numerical evaluation of these eigenvalues has been performed in Matlab using the function eval 
and the parameter values provided for each variant in the main paper, Table 3. 



6.5.2. Finding the eigenvectors. 
The characteristic equation is: 



Av = Xi 



(45) 



or 



-(Pi + G dl + e 12 ) em. -Pi. -Pi 

ei2 -(G d2 + e 2i ) P 2 

G d2 -(P 2 + G r 

Then the eigenvectors satisfying this equation are: 







Vl 




Vl 






v 2 


= Ai 


v 2 


) 




V3 




V3 



(P2+G r +Ai)(Ai + Pi+Gdi+ei2) 
(e 21 -P 1 ){P 2 -G r +\ 1 )-P 1 G d2 

G d2 (Ai+-Pi+G d i+e 12 ) 
(e 2 i-Pi)(P2-G r +A 1 )-P 1 G d2 



1 



(P 2 +G r +A2)(A2+Pi+G dl + e i2) 
(e 2 i-Pi)(P2-G r +A 2 )-P 1 G d2 

G d2 (A 2 +P 1 +G dl +e 12 ) 
(e 2 i-Pi)(P 2 -G r +A 2 )-P!G d2 



(46) 



(P 2 + G,. + A 3 )(A 3 + Pi + G d i+ei 2 ) 
(e 2 i-Pi)(P2-G r + A 3 )-P 1 G d2 

G d2 (A 3 +Pi+G d i+e 12 ) 
(e 2 i-Pi)(P2-G r + A 3 )-P 1 G d2 



(47) 
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6.5.3. The complementary solution. 

The complementary solution can then be written as: 



A 2 t„ 



,A 3 *„ 



y c = C ie Alt v + C 2 e A2t u + C 3 e A3t w 



(48) 



or 



= Ait 



= A 2 i 



=A 3 i 



(P 2 +G r +Ai)(Ai + Pi+G d i+ei 2 ) A] t (P 2 +G r +A 2 )(A 2 +P 1 +G d i+e 12 ) x,t (P 2 +G r + A 3 )(A 3 +P 1 +G d i+ei 2 ) A 3 t 



(eai-POCPa-G.+A^-PiG^ 

G d2 (Ai + P!+G d i+ei 2 ) 
(e 2 i-Pi)(P 2 -G r +A 1 )-P 1 G d2 < 



or, by notation: 



{e 21 -P 1 ){P 2 -G T + \ 2 )-P 1 G d2 

P Ai* G d2 (A 2 +Pi+G d i+ei 2 ) 

(e 2 i-Pi)(P 2 -G r + A 2 )-P 1 G d2 ' 



?y c = 4 C C. 



(e 2 i-Pi)(P 2 -G r + A 3 )-P 1 G d 

= A,t G d2 (A 3 +Pi+G d i+ei 2 ) \ 3 t 

(e 2 i-Pi)(P 2 -G r + A 3 )-P!G d2 e 



c 2 

c, 

(49) 
(50) 



6.5.4- Finding the particular solution 

The following analytical steps are performed in Matlab 
a) We evaluate the inverse of the complementary matrix A 



-l. 

c i 



b) We evaluate the product: 



Z = A7 1 *P = A: 



Pi 





c) We integrate symbolically Z using the commend: R = int(Z,t); 



d) We evaluate the particular solution: y p = A C R. 



(51) 



6.5.5. Finding coefficients Ci, 62,63 by use of initial conditions. 

To coefficients of the homogeneous solution can be found by considering the following condition: 



A c (0)C + y p (0) = y ; 



(52) 
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The above system of equations cannot be solved in Matlab straightforward as the function solve is unable 
to find a solution when all three equations are presented simultaneously; the system can be solved however, 
if the analytical expression of two of the coefficients is given and only one equation is solved symbolically. 
For this purpose, we derive the equations for the first and second coefficient and solve symbolically in Matlab 
only for the last coefficient. We expand the above equation as follows: 



A C11 (0) A C12 (0) A C13 (0) 
A C21 (0) A C22 (0) A C23 (0) 
A C31 (0) A C32 (0) A C33 (0) 





C x 




lfri(0) 









c 2 


+ 


y P2 (o) 


= 







c 3 




y P M 








(53) 



or: 



A C11 (0)C X + A C12 (0)C a + A C13 (0)C 3 + y pi (0) = 
A C21 (0)Ci + A C22 (0)C 2 + A C23 (0)C 3 + y P2 (0) = 
A C31 (0)Ci + A C32 (0)C 2 + A C33 (0)C 3 + y P3 (0) = 



(54) 



We get: 



and 



& = (-y Pl (0)-A Cl2 (0)C 2 -A Cl3 (0)C 3 )/A Cll (0); 



(55) 



C 2 = 



A C21 (0)V P M - A Cii (0)VpM - C 3 (A C23 (0)A Cll (0) - A C21 (0)^ C13 )(0) 



(56) 



A C22 (0)A C11 (0) - A C21 (0)A C12 (0) 

We introduce equations 1551 and [56] in the last equation [54] and then solve symbolically this equation in 
Matlab for C 3 ; we then evaluate C\ and C 2 . 



We are now able to evaluate the full theoretical solution: 



oi = C 1 e Xlt + C 2 e X2t + C 3 e x * t +y pl 

o 2 = A C21 {0)C ie Xlt + A C22 {0)C 2 e X2t + A C23 (0)C 3 e X3t + y P2 ; 



(57) 



The photocurrent elicited in light on condition will be: 
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1 = ^1(01+702). (58) 



or 



I = Vg x [C x {l + 7 A C21 (0))e Al * + C 2 (l + 7 A C22 (0))e A2t + C 3 (l + 7 A C23 (0))e A3t ] + ^ 5l (y pl + 72/p2 ) (59) 
The reader should note that all the eigenvalues above are negative. Therefore we can write: 

Ai = -Ai; A 2 = -A 2 ; A 3 = -A 3 ; (60) 

where Ai, A 2 , A 3 > 0. 

We evaluate the time constants associated with these eigenvalues to be: 

1 l 1 tt-M 

Ti = aT ; T2 = aV T3 = T 3 (61) 

The smallest time constant is associated with the rise phase (from zero to peak) of the photocurrent (the 
the beginning of the optostimulation) ; the other two controls the fast and slow component of the photocurrent 
decay from peak to steady state. With the parameters found in the manuscript Table 3 we identify: 

T rise{on) = T 2 1 T s (on) = T l 5 T f(on) = T 3 (62) 



A rise ( n) = C 2 (l + 7 A C22 (0)); A s(on) = Ci(l + 7 A C21 (0)); A f{on) = C 3 (l + 7 ^ C23 (0)) (63) 



I r ise(on) = V giA rlse ( on y, I s {on) = VgiA s ( on y, If (on) = VgiAf( on y, I p i at = Vgi(y p i + 7 y p2 ) (64) 



Iou{t) = /„ se(on) e- A — c»»)* + / s(on) e- A =(-)* + / /(o „)e- A ^»>' + I p i at ; (65) 



Fo r convenience wc reproduce here the analytical solution for light off condition provided in 
(bfJOJ): 

35 



Nikolic et al 



where 



and 



Ai = 6-c; A 2 = 6 + c; (66) 



Gdl + Gd2 + e-V2 + (321 



b= ^T^T^TC., c = Vb2 ^^^ - Gdie2i - Gd2evi) (6?) 



loff = hiots)^ ' wn + ^/ir /<0//) ; (68) 



^(off) = V0i4,(o//); ^/(o//) = VsiA /(o//) ; (69) 



A s(off) = 7 7 (70) 



[A 2 - (G d i + (1 - 7)ei 2 )]Oio + [(1 - 7)eai + 7(A 2 


- Gd2)]O-20 


A 2 -Ai 

[Gdi + (1 - 7)ei2 - Ai]Oi + [7(Gd2 - Ai) - (1 - 


-7)e 2 i]0 2 o 



A f(°ff) Aa _ Al 
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